from fractions import Fraction
import numpy as np
from sympy import Matrix

#-------------------求行标准形-Hermite矩阵，求满秩分解A=FG，求满秩矩阵的左逆、右逆-------------------
np.set_printoptions(formatter={'all':lambda x: str(Fraction(x).limit_denominator())})#设置np.array用分数表示
#方法：行初等变换[A|I]~[H|P]   (PA=H)
#Herimite的变换矩阵P：对可逆矩阵计算时正确的，不可逆的计算不正确
def hermite_reflection(A):
    rank = np.linalg.matrix_rank(A)
    print("r(A)=", rank)
    I=np.eye(A.shape[0])
    # A_I=np.hstack((A,I))#[A|I]
    # print("A_I:\n",A_I)
    H=np.array(Matrix(np.array(A)).rref()[0].tolist())#A的Hermite矩阵-Matrix().rref()返回两个元素的元组。第一个是精简行梯形形式，第二个是枢轴列索引元组。
    ind_index=np.array(Matrix(np.array(A)).rref()[1])#极大无关组的列索引
    print("极大无关组的列索引:",ind_index)
    # print(A[:,0],A[:,1])
    F=[]
    F=np.array([A[:,i] for i in ind_index ]).T#列极大无关组F
    print("A的满秩分解之列极大无关组F(列满秩矩阵):\n",F)

    G = np.array(H[0:rank, :] ) # 取A的Hermite标准形H的前r行，构成行满秩矩阵G；
    print("A的满秩分解之行满秩矩阵G：\n", G)

    #A=FG
    print("A=FG:",F.dot(G))
    # print("H_I:",H_I)
    print("A的Hermite矩阵H(行标准形):\n",H)
    if rank==A.shape[1] and rank==A.shape[0]:#满秩判断
        A_inv = np.linalg.inv(A)  ##方阵求逆
        # A_inv=np.linalg.pinv(A)##广义逆（也叫“伪逆”）-非方阵求逆  #设置分数显示时报错
        P=np.array(H.dot(A_inv))
        print("A的Hermite变换矩阵P（非行满秩时算的不对）:\n", P)#用广义逆算的不对
    # 列满秩求A的左逆
    if rank==A.shape[1]:
        A_left_inv = np.linalg.pinv(A.T.conjugate().dot(A)).dot(A.T.conjugate())
        print("列满秩矩阵F的左逆：\n", A_left_inv)
    # 行满秩求A的右逆
    if rank == A.shape[0]:
        A_Right_inv = (A.T.conjugate()).dot(np.linalg.pinv(A.dot(A.T.conjugate())))
        print("行满秩矩阵A的右逆：\n", A_Right_inv)
    return H

# A=np.array([[1,2,0,1,5], [2,4,2,8,18], [3,6,1,6,19]])
A = np.array([[1,1,2,4], [-1,-1,-1,-2], [3,3,2,4]])
hermite_reflection(A)
#-------------------求行标准形-Hermite矩阵，求满秩分解A=FG，求满秩矩阵的左逆、右逆-------------------